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ABSTRACT 

We investigate by two-dimensional axisymmetric relativistic hydrodynamical 
simulations (1) jet propagations through an envelope of a rapidly rotating and 
collapsing massive star, which is supposed to be a progenitor of long duration 
gamma ray bursts (GRBs), (2) breakouts and subsequent expansions into stellar 
winds and (3) accompanying photospheric emissions. We find that if the enve- 
lope rotates uniformly almost at the mass shedding limit, its outer part stops 
contracting eventually when the centrifugal force becomes large enough. Then 
another shock wave is formed, propagates outwards and breaks out of the enve- 
lope into the stellar wind. Which breaks out earlier, the jet or the centrifugal 
bounce-induced shock, depends on the timing of jet injection. If the shock break- 
out occurs earlier owing to a later injection, the jet propagation and subsequent 
photospheric emissions are affected substantially. We pay particular attention to 
observational consequences of the difference in the timing of jet injection. We 
calculate optical depths to find the location of photospheres, extracting densi- 
ties and temperatures at appropriate retarded times from the hydrodynamical 
data. We show that the luminosity and observed temperature of the photospheric 
emissions are both much lower than those reported in previous studies. Although 
luminosities are still high enough for GRBs, the observed temperature are lower 
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than the energy at the spectral peak expected by the Yonetoku-relation. This 
may imply that energy exchanges between photons and matter are terminated 
deeper inside or some non-thermal processes are operating to boost photon en- 
ergies. 

Subject headings: black hole physics, hydrodynamics, supernovae, jets and out- 
flows, radiation mechanisms: general 



Introduction 



There is mounting o bservational evidence that links GRBs to the death of massive stars 
(IWoosley fc Blooml 120061 ) and it is widely believed that GRBs are associated wit h the for- 
mation of black hole or magnetar via the collapse of rapidly rotating massive stars (IWoosley 



19931 : |Paczvnskilll998l : MacFadven fc Wooslev 



19991 ). Although we do not know exactly how 



a large amount of energy is generated, the most promising scenario is that a relativistic 
jet is launched from the central engine by neutrino annihilation or magnetohydrodynamical 
proce sses, propagates through a progenitor star and stellar wind (IMacFadyen fc Wooslev 
19991 ). and then dissipates its kinetic energy by internal shocks or pho t ospheric emissions 



or relativistic turbulenc es ( jPiranl 12004 ; iPe'erl 120081 ; lLazzati et al.l 120091 ; lLazar et al.l 12009 
Kumar fc Narayanl 120091 ) . producing the prompt emissions of GRBs or XRFs. 



A large number of numerical works have been d evoted so 



of re l ativistic jet propagations in the stellar envelope (lAlov et al. 



2004 



Mizuta et al. 



Lazzati et al. 



2009 



2006 



Mizuta fc Alov! 120091: Morsonv et al.l 



ar to the understanding 



2000l: Izhang et al 



2007 



Tominaga et al 



2003 
20071 : 



Mizuta et al.l boioh . These simulations have demonstrated that the jet 



is confined by the pressure of hot cocoon as it penetrates through the stellar envelope. 
The Kelvin-Helmholtz in stability, which occurs between the cocoon and the jet, produces 



rich internal s tructu r es (lLazzati fc Bege 



man 



20051 : Morsonv et al.l 120071 ). More recently 



Lazzati et al.l ( 120091 ): iMizuta et al.l ( 120101 ) computed the jet propagation beyond the stel- 
lar surface and observed that these internal structures in the jet and cocoon leave their 
traces until later times. They also claimed that the hot jet produces very bright and highly 
efficient photospheric emissions in the prompt phase of GRBs. These very effi cient pho- 



tosph eric emissions may solve the efficiency problem of the prompt emission ( lloka et al. 



2007 1- Interestingly, thermal emissions were indeed ide ntified for some Long GRBs lately 



(Abdo et al. 


2009; 


Rvde et al. 


2010; 


Guiriec et al. 


2010) 



2008; 


Toma et al. 


2010; 


Ioka 


2010) 
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It should be noted that the previous numerical studies on the jet propagation ignored 
the infall of the stellar envelope. According to the collapsar model, on which this paper is 
based, the gravitational core collapse sets in just like ordinary supernovae when the density 
reaches p ~ 10 9,5 g/cm 3 or the temperature exceeds T ~ 5 x 10 9 K and electron captur es or 



endothermic photodissociations of nuclei reduce pressure (see e.g. iKotake et al.ll2006l ). A 
shock wave produced by core bounce stalls in the core and a large amount of matter accretes 
on a time scale of seconds onto a proto-neutron star at first and into a black hole later. The 
so-called prompt shock wave either remains stagnated near the black hole or is swallowed 
into it. On the other hand, the core collapse produces a rarefaction wave at the boundary of 
the core and envelope, which then propagates outward through the envelope and induces the 
infall of the envelope when it arrives. Thus, the neglect of the envelope motion in studying 
the jet propagation is justified only when the jet is launched very early o n, possibly soon 



after the black hole is formed (IMacFadyen et al.l l200ll ; IZhang et al.l 120031 ) , and the infall 



of the envelope is not yet substantial. If the jet launch is delayed somehow, on the other 
hand, the profile of the envelope will be modified and the jet propagation will be affected. 
It is also pointed out that the stellar envelope may cease to infall eventually. In fact, the 
outer portion of the stellar envelope is likel y to have an angular m omentu m large enough 
to ter minate the infall by centrifugal forces (jWoosley Sz Hege 3 hood ). Indeed, iLindner et al.l 



(120101 ) observed in their long term simulations of rotational collapse of massive stars that 
a shock wave is generated by centrifugal forces and the outer portion of stellar envelope is 
expelled eventually. Since we do not know exactly when the jet is launched, it is important 
to study the influence of envelope dynamics on the jet propagation and subsequent prompt 
emissions. 

Motivated by these facts, we numerically investigate the relativistic jet propagations 
through a non-stationary envelope, moving either inward or outward, of a rapidly rotating 
massive star, varying the timing of jet injection. We assume that the prompt shock wave 
of core bounce origin has already been swallowed into the black hole and what is supposed 
to occur in the core after bounce such as neutrino heating and various hydrodynamical 
instabilities do not affect the dynamics of the envelope. We do not specify the mechanism of 
the jet launch from the central engine, which is under controversy at present, and inject the 
jet with appropriate properties from the computational inner boundary by hand, following 
the common practice in this field. Our focus is the jet propagations in the non-stationary 
envelope and its influences on the subsequent photospheric emissions. This paper is organized 
as follows. In Section 2, we describe the models and numerical procedures. Then our main 
results will be presented in Section 3. We conclude the paper with the summary of our 
findings in Section 4. 
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2. Methods 

As mentioned above, in this letter we compute the jet propagations through the envelope 
of a rotating massive star into a stellar wind, taking into account the core-collapse-induced 
motions of the envelope under the assumption that the prompt shock wave is soon sucked 
into the black hole and various processes in the core such as neutrino heating of accreting 
matter and hydrodynamical instabilities do not affect the dynamics of the envelope. In order 
to simulate the infall of the envelope induced by the rarefaction wave that is generated by 
core collapse, we take rather involved multiple steps. More specifically, we (1) construct 
massive star's envelope models in rotational equilibrium, (2) put quasi-steady winds on top 
of them, (3) simulate rotational collapse of the envelope, generating a rarefaction wave by 
artificially reducing the pressure gradient at the inner boundary and (4) compute subsequent 
jet propagations in the envelope and (5) calculate photospheric emissions as a postprocess. 
We employ the so-called HSCF scheme in the first step and perform 2D relativistic hydro- 
dynamical simulations in the third and fourth steps. In the following, we explain what is 
done at each step more in detail in order to facilitate reader's understanding of our results 
in the next section. 



2.1. A Massive Star Envelope in Rotational Equilibrium 

The first step is a preparation of the initial model for dynamical simulations in the 
later phases. In this subsection, we construct a 2D axisymmetric model of a rotat ing mas- 



(1986 


); 


Kiuchi et al. 


(2010) 



evolution models are still unable to fully implement rotational equilibrium and neglect the 
non-spherical deformation of rotating stars. In this study, however, the rotational equilib- 
rium is crucial, since the infall of the envelope commences only after the rarefaction wave 
generated at the boundary of the core and envelope arrives. If the initial model is not in 
dynamical equilibrium, however, even outer parts of the envelope begin to move immediately 
after the simulation is started and false shock waves are produced as a consequence more 
often than not. 



Our envelope model is constructed so as to mimic 16TI model by IWoosley fc Heger 



( 120061 ). which is currently supposed to be one of the most promising GRB progenitor models. 
Since the outer envelope of 16TI is almost radiation dominated, we employ a polytropic 
equation of state (EOS) with the adiabatic index of 7 = 4/3. We impose a rigid rotation as 
an approximation to the outer envelope of 16TI. Figure [T] shows the density profiles on the 
rotational axis and equator for our model together with the one for 16TI. Also displayed in 
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the figure is the density distribution in the meridian section for our model. Our model agrees 
fairly well with 16TI except for the innermost portion, where the rotation is not rigid in 16TI. 
This discrepancy is not very important for the investigation in this study, since that part 
is sucked into the inner boundary much earlier than the jet injection. Our envelope model 
has a total mass of M ~ 14M Q and a specific angular momentum of j sp ~ 1.5 X 10 19 cm 2 /s 
at the stellar surface, which is close to the mass shedding limit. The rotation velocities of 
our model are slightly lower than those of 16TI in general. The specific angular momentum 
distribution on the equatorial plane as a function of enclosed mass is shown in Figure [2J 
Here the enclosed mass is defined as a mass within a certain radius. Presented also in this 
figure is the specific angular momenta at the innermost stable circular orbit (ISCO) for 
Schwartzshicld black holes as a function of their masses. The two curves intersects with each 
other at the enclosed mass of ~ 8M . Envelope matter that has a larger enclosed mass than 
this value can not fall down to the black hole and is halted somewhere outside the ISCO by 
centrifugal forces. As a matter of fact, we find a centrifugal bounce and a formation and 
subsequent propagation of a shock wave (see subsection 12.41) . 



2.2. Special Relativistic Hydrodynamic Code 

Before proceeding to the subsequent steps, we describe the numerical methods employed 
for the hydrodynamic simulations done in Steps 2 to 4. We employ a 2D axisymmetric, special 
relativistic hydrodynamics code. Equatorial symmetry is also assumed in this paper. The 
basic equations we solve in this study are given as follows in the geometrical units G = c = 1, 
where G and c are the gravitational constant and speed of light: 

dtp* + dj (p*v J ) = 0, (1) 
d t S r + dj (r 2 sin 6 T\) = r 2 sin #{-T 00 W + r T ee + r sin 2 9 T^j, (2) 

d t S + dj(r 2 sin9 T j ) = r 2 sin ^j-T 00 ^ + r 2 sin 9 cos 6T**\ (3) 

dtSt + dj^wxeT*^ =0, (4) 
d t r + dj (r 2 sin 9 T 0j - pj?) = -r 2 sin 9 T°V,i, (5) 
d t (p*A) + dj (p*Avi) = 0, (6) 

= 47tp hh(u t ) 2 -h + 2^-\, (7) 

where the subscript j runs over r and 9, and A, T^, and ip denote the mean molecular 
weight, energy-momentum tensor of ideal fluid, four- velocity of matter and gravitational 
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potential, respectively, and 

= r 2 sinOpou 1 , (8) 

Si = r 2 smeTl (9) 

r = r 2 sin 6T 00 - p*. (10) 

The above equations are derived from the Einstein equations and energy-momentum con- 
servation equations by the weak field approximation, ignoring the time-derivative of gravi- 
tational potential and space-derivatives of three dimensional space metric. Since our com- 
putational domains do not contain the origin, the gravity of the central object is added as 
a point mass at the center. The time evolution of the mass of the central object is taken 
into account by integrating the mass flux crossing the inner boundary of the computational 
domain. 

We solve the Poisson equation for the gravitational potential, Eq . (171), by MICCG and 



the hydrodynamical equations, E qs. flTJ-flB]), by the central scheme (jKurganov fc Tadmor 



20001 ; iNagakura fc Yamadal 120081 ) . In the latter, using the PPM interpolation method and 
TVD Runge-Kutta time integration, we achieve the second order accuracy in both space and 
time. 

The EOS's employed in this paper are the following. For Steps 2 and 3, that is, the 
construction of t he ste llar wind and the computation of the envelope collapse, the EOS by 



Blinnikov et al.l ( 119961 ) is used, in which the temperature and mean molecular weight are 
introduced to avoid the inconsistency with Step 1, where they are also accounted for. On 
the other hand, the so-called 7- law EOS, p = (7 — l)poe with p, 7 = 4/3, po and e being 
the pressure, adiabatic index, rest mass density and specific internal energy, respectively, is 
adopted for the jet simulations in Step 4 for simplicity. Since we find that the envelope is 
radiation-dominated at the time of jet launch, this is a good approximation. 

It is a consensus that high resolution simulations are necessary for the investigation of 
interactions between the jet and stellar envelope in the jet drilling phase, since the Kelvin- 
Helmholtz instability and turbulent motions take place inevitably. When the velocity of jet 
head is smaller than the local sound speed at the hot spot, which is indeed the case for 
the jet propagation s in the stellar envelope, a back flow is bent and pinches the jet path 



( iMizuta et al.ll2010l ). In order to treat these effects adequately, we employ an adaptive mesh 
refinement (AMR) technique, in which the forward shock is searched at each time step and 
the number of mesh points in its vicinity is increased in each coordinate direction. 

In Appendix, we show the results of several numerical tests meant to validate our hydro- 
dynamics code used in this paper. It is also demonstrated that the rotational massive stellar 
envelope, which is constructed by HSCF scheme at Stepl, does not change the configurations 
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in a dynamical simulation, which is clear evidence the both codes are reliable. 



2.3. A Quasi-Steady Wind 

Massiv e stars experience ma s s losses in general and the GRB progenitors will not be 



exceptions ( jCampana et al.ll2006l ; ISoderberg et al.l 120081 ) . We hence take into account the 



stellar wind in our initial model (Step2). It is noted, however, that theoretical understand- 
ing of the driving mechanism of stellar winds and mass losses of massive stars is far from 
satisfactory and it is much beyond the scope of this paper to address these issues. We are, 
therefore, satisfied with the construction of quasi-steady winds, not specifying its driving 
mechanism. It is stressed that what is important here is that the wind thus obtained does 
not change its configuration very much before the jet reaches it. 

We first construct a spherically symmetric, steady wind configuration, neglecting rota- 
tion. For this purpose, we perform ID hydrodynamical simulation, using the code described 
above, in the region from the stellar surface up to the distance of 10 13 cm. The initial config- 
uration is rather arbitrary. Fixing the density, pressure, velocity (or mass loss rate) at the 
inner boundary, we run a long-term simulation until the wind is settled to a steady state. 
The values of the density and pressure are chosen so that they would be continuous when 
the wind is appended to the envelope model constructed at Step 1. Rotation is then added 
so that the specific angular momentum should be constant along each radial ray. The values 
of the specific angular momentum at the inner boundary are chosen in such a way that they 
would be continuous from the envelope to the wind. The wind obtained in this way is not 
exactly steady any more. Although rotational, steady wind configurations could be obtained 
in the similar way, it turns out that the wind configuration does not change much during 
the jet propagation through the envelope and wind if rotation is added this way. We hence 
do not pursue further elaboration in this paper. 

By changing the inner boundary condition, we can construct various wind models, both 
dense and tenuous. In this paper, however, we adopt only an optically thin model to elucidate 
the effects of envelope motions on the jet dynamics. Other wind mod els and their influence s 



on the jet propagation will be investigated in the sequel to this paper (INagakura et al.ll201ll ). 
The photosphere of the present wind model is located at the stellar surface and the mass loss 
rate is M ~ 10~ 6 M Q /yr. Figure [3] shows the profiles of our wind model. The density and 
pressure distributions obey power laws approximately with the power- law indices being —2.14 
and —2.82, respectively. The outflow in the wind becomes supersonic at r ~ 7.5 x 10 10 cm 
and its velocity approaches asymptotically v r asym ~2x 10 8 cm/s. 
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2.4. Collapse of the Massive Star Envelope 

Using the envelope and wind configurations obtained above as an initial condition, we 
perform 2D axisymmetric simulations of the envelope collapse (Step3). The computational 
domain covers at first a region of 5 x 10 6 cm < r < 2 x 10 12 cm, which includes the entire 
envelope and the core region except the black hole and its close vicinity as well as the inner 
part of the wind. The inner boundary is shifted outwards later (see below). The radial grid 
consisting of 224 points is non-uniform, with the grid width changing with the density scale 
height. The angular grid covers a quadrant of the meridian section and is uniform with 60 
points. 

In reality, as we mentioned earlier, the gravitational collapse of the envelope is initiated 
by the arrival of the rarefaction wave that is generated at the core/envelope boundary by 
core collapse and propagates outward. To mimic this, we reduce the radial gradients of all 
quantities to zero at the inner boundary and artificially induce the infall there. Then a 
rarefaction wave is produced at the inner boundary and propagates outward, inducing infall 
at points it reaches. It is stressed that we confirmed by the test computation presented 
in Appendix IA.6I that if we do not reduce the radial gradients of quantities at the inner 
boundary, the envelope remains intact even after many time steps. 

As shown in the next section, the contraction of the envelope is eventually terminated by 
centrifugal forces, producing a shock wave that propagates outwards and eventually breaks 
out of the stellar envelope. We increase the number of radial grid points to 1000 at the time 
of the shock breakout and shift the inner boundary outwards to 5 x 10 8 cm simultaneously. 
All the quantities are linearly interpolated to the new mesh points. The change of the inner 
boundary leads to the increase of the mass of the central object, which is properly taken into 
account, whereas we discard the angular momentum and energy between the old and new 
inner boundaries just for simplicity. 

It should be noted that our numerical code does not take into account general relativity 
and detailed microphysics such as photo-dissociations of nuclei and neutrino cooling. The 
neglect of these effects tends to overestimate the strength of the shock wave of centrifuga l 



bounce-origin. In fact, it was pointed out by ( ILindner et al.ll2010t iMilosavljevic et al.ll2010f ) 



that the nuclear photo-dissociations may completely sap the shock wave. We will defer the 
investigation of this issue to a future work, in which we will implement a nuclear network in 
our hydordynamics code. It is also repeated that we assume in this paper that the prompt 
shock wave of core bounce-origin is swallowed into the central black hole and what occurs 
inside it does not affect the dynamics of the envelope. In order to see if this assumption is 
correct or not, it is necessary to perform detailed simulations of core collapse in full general 
relativity, which is a major undertaking and will also be a future work. 
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2.5. Jet Injection and Propagations through the Stellar Envelope and Wind 

In the next step (Step4), which is the main part in this paper, we numerically study 
the jet propagations through the stellar envelope and wind that are in motion as obtained in 
the previous step. Following the common practice, we inject a relativistic jet from the inner 
boundary, not specifying the driving mechanism, at two different times after the envelope 
collapse takes place: 20s for model M20s and 50s for model M50s. The injection parameters 
are identical for both models: the jet is hot (p/poc 2 = 20, where c is the speed of light) 
and relativistic with a Lorentz factor of 5; the half opening angle is 9°; the power of jet is 
constant in time and the injection continues for t<i ur = 30s with the total injected energy 
being 10 53 ergs. Then the terminal Lorentz factor is estimated by 

Tfer-m = h in T in = (1 + € in + Pin/ (PinC 2 ))T in ~ 7/(7 - l) X p in /(p in C 2 ) X T in , 

= 4{Pin/(PinC 2 )}T in , (11) 

where h in , p in , p in and Y in are the specific enthalpy, pressure, rest mass density and Lorentz 
factor at the injection; the adiabatic index is denoted by 7 and is set to be 4/3. The choice 
of the injection parameters in this paper corresponds to T term ~ 400. 

The computational domain for these simulations ranges from r = 10 9 cm up tp r = 
10 18 cm. Note that this broad range is mandatory for the identification of the location of 
photosphere until t a 5 S ~ 100s, since the forward shock in the jet is highly relativistic with 
a Lorenz factor of T > 100. The total number of radial grid points is 11000. The grid is 
nonuniform, with the grid width being smallest (Ar = 10 8 cm) at the inner boundary and 
increasing geometrically by ~ 0.1% per zone up to 10 13 cm and by ~ 1.35% in the region 
further out. The number of angular grid points is the same, 60, as in the previous step. We 
remap the data obtained in the previous step to the new grid by the same linear interpolation 
as employed in Step3. The shift of the inner boundary requires an adjustment of the mass 
of central objects, with the mass between 5 x 10 6 (5 x 10 8 )cm < r < 10 9 cm being added to 
the central point mass for model M20s (M50s). The density, pressure and velocity in the 
region of 10 13 cm to 10 18 cm are extrapolated from the inner region in the following manner: 
The density and pressure are extended by the power-laws that fit their distributions in the 
inner region; The radial velocity is assumed to be constant in the extended region, since 
it has already reached the asymptotic velocity (see the bottom panel of Figure [3]); The 9 
component of velocity is set to be 0, whereas the azimuthal component is determined so that 
the specific angular momentum is constant along each radial ray just as in Step 2. 

During the jet propagation through the stellar envelope, we employ an AMR technique. 
In our code there are only two levels of meshes deployed, in which the resolution of the 
second level can be varied. Here the mesh of the second level is 9 times finer than the first 
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level mesh, with the smallest radial and angular resolutions being Ar = 1.1 x 10 7 cm and 
A8 = 0.16°, respectively. After the jet breakout, on the other hand, the jet head expands 
nearly freely and soon becomes highly relativistic. As a result, the back flow tends to be 
suppressed and the jet morphology does not change so much during this phase. We hence 
employ only 3 times finer a mesh for the second level after the jet head reaches R = 10 11 cm . 
The resolution in this study is not as high as in the previous study (ILazzati et al.l 120091 ). 
One of the main reasons for this is the fact that we are dealing with a much greater spatial 
extent. This is necessary, as already mentioned, to identify the locations of photosphere. As 
a result, however, rapid variations in the photospheric emissions are sacrificed to some extent 
by numerical dissipations and our discussions on this issue are restricted to a qualitative level. 



2.6. Photospheric Emissions 

As a final step (Step5), we calculate, as a post process, the photospheric emissions based 
on the data obtained in Step4. We define the photosphere to be the surface that has a unit 
optical depth from infinity with respect to the Thomson scattering. The optical depth is 
given by 

/■oo 

r(t obs ,r)= / n e (t\s)a T T(t*,s)(l- P(t*,s)cos6 v (t*,s))ds, (12) 

J r 

where s is the distance along the line of sight, n e is the number density of electron in the 
comoving frame, ctt is the cross section of the Thompson scattering, /3 is the matter velocity 
normalized by the speed of light c, T is the correspondi ng bulk Lorentz factor an d 9 V is the 



angle between the line of sight and the matter velocity (lAbramowicz et al.lll99ll ). It should 
be stressed that the time retardation expressed by t* = t a b s — s/c in the above equation cannot 
be ignored for relativistic flows. In this paper, we evaluate Eq. (|T2l as it is, retrieving the 
data for appropriate times from the results of the hydrodynamical simulations. Thanks to 
the wide spatial range of our hydrodynamic simulations, photons observed at t a b s < 100s 
have passed the forward shock by the end of simulations, the fact which is important for the 
identification of the locations of photosphere. 

The observed isotropic luminosity of photospheric emissions is then given by 

L = 4vr J V 4 1 cos 6 ph dS. (13) 

Here dS is the areal element of the photosphere (measured in the laboratory frame), T> = 
[r(l — ,5 cos 6*„)] _1 is the Doppler factor, # p h is the angle between the line of sight and 
the normal vector of the photosphere, / = <7sbT 4 /7t is the radiation intensity with a^B 
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and T being the Stefan-Boltzmann's constant and the temperature in the comoving frame, 
respectively. We ignore the cosmological redshift in this study. 



Results 



3.1. Envelope Collapse 



The upper panel of Figure H] shows the temporal evolution of the density profile on the 
equator obtained in Step3, that is, the computation of envelope collapse. The infall starts at 
the inner boundary, generating a rarefaction wave that propagates outwards. Only after this 
rarefaction wave arrives, other parts of the envelope begin to move inwards. The contraction 
is almost spherical initially. As time passes and more distant portions of the envelope start 
to infall, however, the centrifugal force becomes non-negligible, since the specific angular 
momentum is an increasing function of radius. The centrifugal force becomes large enough 
eventually to halt the infall of matter and a shock wave is generated. This happens at 
t ~ 18 s in our m o del. A similar but a bit later bounce by centrifugal force was also reported 
by iLindner et al.l ( 120101 ) . The reason why we found the earlier formation of the shoc k wave is 
that we put the inner boundary at much smaller a radius than lLindner et al.l (120101 ) . Indeed, 
the inner boundary of our model is initially l ocated at ~ 3 times the Schwarzschild radius, 



which is 10 times smaller than that adopted in lLindner et al.l (l2010f ). It should be noted that, 
as we have already mentioned, the shock wave is expected to be produced near the innermost 
stable circular orbit in reality and more accurate computation of its formation requires an 
implementation of general relativity as well as micophysics such as neutrino transports and 
photo-disintegrations of nuclei, which may sap the shock wave. 

The shock wave propagates more vigorously along the equator than along the rotational 
axis and reaches the stellar surface on the equator at t ~ 31 s (see the bottom left panel of 
Figure H|). Then the shock wave breaks out of the stellar surface and runs further through 
the wind (see the bottom right panel of Figure H]). If the jet is launched earlier than the 
shock formation, the shock dynamics just described will be modified by the jet propagation. 
If the opposite is true, that is, the jet launch is later than the shock formation, the jet 
dynamics will be affected by the shock propagation. In particular, if the jet launch is suffi- 
ciently delayed, the jet propagation in the wind and, as a result, the photospheric emissions 
will be severely changed. Model M20s is meant for the former case whereas the latter case 
corresponds to Model M50s. Incidentally, the shock breakout in the latter case may account 
for the so-called precursor that is observed for some long GRBs. In fact, t he typical tim e 
lag between the precursor and the prompt emission is several tens of seconds (jLazzatill2005l ). 
which is similar to what we find in our model. In this scenario, the high energy emissions 
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in the pre cursor are supposed to be similar t o those in the shock breakout of ordinary su- 



per novae (Falk 



Soderberg et al 



978 



Klein fc Chevalier 1978; Matzner fe McKee 1999; Waxman et al. 2007 



20081 ). For more quantitati ve arguments for the precursor emissions from 



these shock waves are currently undertaken (INagakura et al.ll201ll ) 



Incidentally, we assume in this paper that it is not the centrifugal bounce-originated 
shock wave but something else that is responsible for the jet launch. We hence treat the 
centrifugal bounce and the jet launch as independent events and vary the time of jet injection 
with respective to the centrifugal bounce rather freely. In reality, they may be correlated 
with each other one way or another. As we have already mentioned, the focus in this paper is 
the consequences that the possible time lag between these two events may have. The origin 
of the lag is intimately related with the mechanism of the jet launch. Although it is very 
interesting in its own right, the issue is much beyond the scope of this paper. 



3.2. Jet propagations in the Stellar Envelope and Wind 

As expected, the hydrodynamics of the early injection model, M20s, (see the left column 



of Fi gure [5]) is similar to those found in the previous studies (ILazzati et al 



2009; Mizuta et al. 



20101 ). since the envelope bounce by centrifugal force occurs almost at the same time as the 
jet launch and the outer envelope structure has not been changed very much from the initial 
one. The jet is strongly collimated by a hot cocoon, i.e. the shocked jet and envelope 
matter, until the jet breaks out of the progenitor surface. Then the shocked jet matter starts 
to expand laterally from the vicinity of rotation axis and the internal energy is gradually 
converted to the kinetic energy. As a result, the hot, shocked jet matter acquires a high 
Lorentz factor and produces very bright photospheric emissions (see next subsection). Since 
the jet injection is terminated at t = 30s in this model, a rarefaction wave is generated at 
that point and starts to chase the jet head and only the matter between the jet head and 
rarefaction wave contributes to subsequent radiations. 

For model M50s, in which the jet is launched much later than the envelope bounce, the 
jet dynamics is very different from the one for the early injection case (see the right column 
of Figure [5]). The jet propagates through the envelope that is not contracting but expanding 
owing to the shock wave produced at the centrifugal bounce of envelope. We find that the 
distance between the terminal (reverse) and forward shocks is shorter than for model M20s 
and the terminal shock remains to exist much longer for model M50s. The forward shock 
region in the jet is also found to be remarkably different after the breakout between the 
two models. Since the shock wave breaks out of the star before the relativistic jet reaches 
the stellar surface, the stellar wind is modified substantially by the shocked envelope matter 
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(SEM). As a consequence, the jet propagation is hindered by the thick SEM even after its 
passing the position of the original stellar surface and the forward shock velocity becomes 
slower until much later times when the jet passes all through the SEM, producing a denser 
shell behind the forward shock (see the second panel in the right column of Figure E]). This 
has an important ramification for the photospheric emissions later on. 



3.3. Photospheric Emissions 

In Figure E] we display the light curves together with the evolutions of photospheric radii 
and observed temperatures (T>T) for both models. The observer is assumed to be located on 
the rotational axis (on-axis observer) and the photospheric radius in the figure is the value 
on the axis although we calculate the positions of photosphere for off-axis rays. For model 
M20s, the luminosity is peaked in an early phase (t b s < 10s) and rather high luminosities 
are sustained for the next 30s whereas a strong peak is observed at a late time (t Q b s ~ 25s) 
for model M50s. As already mentioned, this difference arises from the large difference in 
the envelope structures prior to the jet breakout. Since a larger amount of matter is swept 
up by the jet in model M50s, it takes the photosphere longer to leave the forward shock 
region and move inward to a region with higher observed temperatures, producing bright 
radiations. This qualitative difference in the light curves may be utilized to observationally 
extract information on the timing of jet launch at the central engine. 

As an explanation of the GRB prompt emissions, the photospheric emissions in our 
models have sufficiently high luminosities. The peak energy, however, is l ower roughly by 



an or der of magnitude than the value expected from the Yonetoku-relation (jYonetoku et al. 



20041 ). This tendency is the main drawback of the photospheric emission model. Note, 
however, that the shocked jet matter may be scattering-dominant and energy exchanges 
between photons and matter may be terminated deeper inside. If this is the case, the observed 
tempe ratures will be h igher and the luminosity will be also reduced. They are currently under 



study (llto et al.ll201ll ). It is also conceivable that some non-thermal processes are operating 
to produce high energy photons. Further exploration of these issues will require detailed 
computations of radiation transport and will also be a future work. 



4. Summary 



We have numerically investigated the propagations through a rapidly rotating massive 
star envelope of relativistic jets that are launched at different times, taking into account the 
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motions of the envelope induced by core collapse. Then, we have calculated the photospheric 
emissions by post-processing. The main findings in this paper are summarized as follows: 

1. In the envelope collapse, we have seen the generation of a shock wave by centrifugal 
force around t ~ 20s for the progenitor rotating uniformly at the mass shedding limit. In 
~ 10s the shock wave breaks out of the star if the jet launch is sufficiently delayed. 

2. If the shock wave produced by centrifugal force breaks out of the star earlier than 
the jet does, it changes the envelope and wind structures drastically and the jet propagation 
thereafter is also affected significantly. In fact, since the forward shock in the jet sweeps up 
a larger amount of matter, a dense shell is produced behind it in that case. 

3. The light curve of photospheric emissions is qualitatively different if the jet is launched 
later and propagates in the shock-modified envelope and wind. In the case of earlier launch, 
the peak luminosity is attained at a relatively early time (t b s ~ 10s), whereas it takes longer 
(tabs ~ 25s) to observe the peak for the delayed launch case owing to the dense shell just 
mentioned. 

4. The photospheric emissions obtained in this study with the time retardation be- 
ing taken into account appropriately have high luminosities suitable for the GRB prompt 
emissions. However, the peak energy tends to be lower than expected from the Yonetoku- 
relation. If the shocked jet matter is scattering-dominated, photons will cease to exchange 
energy with matter deeper inside, where the temperature is higher. It is also possible that 
some non-thermal processes boost the photon energy. 



Lazzati et al. 


( 


2009 


); 


Mizuta et al. 


( 


2010) 



possible consequences that the difference in the timing of jet launch may have for the prompt 
emissions, we have not made detailed comparisons with these previous studies. There are 
a couple of conceivable causes for the differences: (1) better estimation of the location and 
temperature of photosphere thanks to the wider computational domain, (2) the effect of 
envelope collapse taken into account, (3) the differences in the jet-injection parameters, pro- 
genitor models and numerical resolutions. These issues will be addressed in our forthcoming 
paper. 



We thank Akira Mizuta and Shigehiro Nagataki for useful discussions. This work is par- 
tially supported by the Japan Society for Promotion of Science (JSPS) Research Fellowships, 
Grant-in-Aid for the Scientific Research from the Ministry of Education, Culture, Sports, 
Science and Technology, Japan [Nos. 222913, 22740178, 21540281, 19104006]. This study is 
also supported by Program for Improvement of Research Environment for Young Researchers 



- 15 - 



from Special Coordination Funds for Promoting Science and Technology (SCF) commissioned 
by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. 



A. Code Tests 

In this appendix, we carry out a series of tests in order to validate our special relativis- 
tic hydrodynamics code, which employs the PPM reconstruction and TVD Runge-Kutta 
integration with a second order accuracy in both space and time. We adopt a HLL-type 
numerical flux and the CFL number is set to be 0.5. Included in the following are (1) ID 
special relativistic shock tube problems, (2) the same as (1) but with tangential velocities 
and (3) a 2D Riemann problem. The ID problems are compared with the exact solutions, 
and we utilize the results given in previous papers for the 2D problem. We also solve (4) 
one-dimensional and (5) two-dimensional isentropic flows to obtain the convergence rate 
quantitatively. For these test runs (l)-(5), the 7-law EOS is adopted with the adiabatic 
index of 7 = 5/3. (6) We also run a dynamical simulation of the rotating stellar envelope, 
which is obtained by the HSCF method (see section 12.11) and confirm that the stellar enve- 
lope sustains the initial profile for a long time. In order to check the accuracy of our AMR 
part, (7) we compute a non-relativistic, spherical point explosion, which can be compared 
with the Sedov- Taylor analytical solution, (8) a pulse propagating adiabatically through 
meshes of different refinement levels, and (9) an axisymmetric, relativistic jet propagation 
in a uniform matter. These tests demonstrate that our numerical code has enough accuracy 
for the purpose of the current study. Throughout this appendix, we adopt geometrical units 
G = c — 1 otherwise stated. 



A.l. ID Relativistic Shock Tube Problems Without Tangential Velocities 

The shock tube problem is one of the common tests for hydrodynamical codes. It is 
a special Riemann problem in gas dynamics. One of t he advantages of this test is the fact 



that we know exact solutions even in special relativity ( jPons et al.ll2000l ). We can check how 
well the code reproduces the profile of a rarefaction wave and captures several discontinuities 
such as contact surface and shock wave. In this test, we set the number of grid points to be 
400 and the parameters employed for two runs are as follows: 

case 1. Left state: (p,v,p) L = (10,0,13.3), Right state: (p,v,p) R = (1,0, 10" 6 ) 

case 2. Left state: (p,v,p) L = (1,0, 10 3 ), Right state: (p,v,p) R = (1,0, 10" 2 ) 

Figure [7] shows the results at t = 0.4 together with the exact solutions. As is obvi- 
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ous, the overall profiles are well reproduced. Although the contact surface and shock wave 
are somewhat smeared out, ou r results are quite similar to those of other groups (see e.g. 



Del Zanna fc Bucciantinill2002l ). 



A. 2. ID Relativistic Shock Tube Problems With Tangential Velocities 



Here we show the results of relativistic shock tube problems with tangential velocities. In 
the special relativistic shock tube problems, the velocity components tangential to a discon- 
tinuity play a non-trivial role unlike in the non-relativistic counter part because the Lorentz 
factor depends on the absolute value of the velocity a nd it is numerica l ly harder to resolve 



the f low profiles in special relativity as reported by (IPons et al.l |2000| ; iRezzolla fc Zanotti 



20021 ). In these tests, we adopt the same initial condition as in case 2 of the previous sub- 
section excep t for t he non-vanishing tangential velocities, which are identical to those in 
Mizuta et all J2006h . 



Figure [8] shows the results of these tests. We have varied the tangential velocity v y 
from to 0.99 on both sides. It is clear that both the contact surface and shock wave are 
substantially deviated from the exact solutions as the tangential velocity becomes large. We 



have performed test runs for (v^. 



(0.9, 0.9) with higher spatial resolutions (the number 



of grid points change from 800 to 6400) and display the results in Figure EO Although the 
deviations of the numerical results from the exact solution are still noticeable even in these 
high resolution runs, we can confirm the convergence of the numerical results to the exact 



solution. Again the perform ance of our code is similar to others (see e.g. iMizuta et al.l 12006 
Zhang fc MacFadverJ [20061 ) . 



A. 3. A 2D Riemann Problem 

This test is meant to check the performance of our code in multi-dimensional set- 
tings. The computational domain is initially divided in 4 sections, which have different 
states. The solution consists of multiple shock waves, contact surfaces and a rarefaction 
wave interacting with each othe r. The parame t ers w e adopt in this test are the same as in 
Del Zanna fc Bucciantini] (120021 ): IMizuta et al.l kood \: iMorsonv et aD ((2003): 



(p,v x ,v y ,p) = (0.1, 0.00,, 0.00, 0.01) for 0.5 > x > 1, 0.5 > y > 1 
{p,v x ,v y ,p) = (0.1,0.99, ,0.00, 1.00) for > x > 0.5, 0.5 > y > 1 
{p,v x ,v y ,p) = (0.5,0.00, ,0.00, 1.00) for > x > 0.5, > y > 0.5 
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(p,v x ,v y ,p) = (0.1,0.00, ,0.99, 1.00) for 0.5 > x > 1, >y > 0.5 

We use a uniform mesh with 400 x 400 grid points. Figure [TD] shows a contour in 
the logarithm scale of the rest mass density obtained in this simulation. There is no exact 
solution available. The results appear very similar to the ones presented in previous studies. 



A. 4. A ID Isentropic Flow 

All the above tests involve discontinuities such as shock wave and, as a result, the 
code affords only a first order accuracy because the numerical error is dominated by these 
structures. Note that this is necessary to ensure numerical robustness. In order to see the 
performance of our code for smooth flows, we have carried out numerical simulations of ID 
and 2D (see the next subsection) isentropic flows. The exact solutions are obtained by the 
characteristic method. The test hence offers an opportunity to quantitatively assess the 
accuracy of our code. 

The initial condition for this test is the same as those employed in the previous studies 



(jZhang &: MacFadyenl 120061 ; iMorsony et all 120071 ) and the density profile is given by 



Po(x) = pref{l + af(x)}, (Al) 
where p re f is the density of a reference state and 

nX) -\0 : otherwise (A2) 

with a and L being the amplitude and width of a pulse, respectively. Since the flow is 
isentropic, we employ a polytropic EOS (p = Kpl) with the polytropic constant of K — 100 
and the adiabatic index of 7 = 5/3. The velocity of the reference state is set to be 0, while 
the velocity distribution inside the pulse is chosen so that the left-going Riemann invariant 
should be constant. Thanks to this set-up, the wave propagates in one direction. The special 
relativistic Riemann invariants are given by 



J ± = ±\n(—) ± -^M^==^), (A3) 

where c s denotes the sound velocity. The equations of characteristics C± are expressed as 

dx\ v±c s (A4) 



dt J G 1 ± vc s 
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Since the J_ and entropy are constant over the whole region, the pulse evolution is determined 
by the J + , which is carried along the characteristic C + until right-traveling characteristics 
collide with each other and a shock wave forms. Although the shape of pulse is initially 
symmetric, it is skewed owing to different characteristic velocities. Note also that the post- 
pulse state is the reference state. 



Our computational d omain is the same as in previous studies (IZhang fc MacFadyen 



20061 ; iMorsony et all 120071 ): (—0.35 < x < 1); The reference state has p re f = 1, p re f = 100 
and v re f = 0; The amplitude of the pulse is a = 1.0 and the width is L = 0.3. The 
simulation is run until t = 0.8. A comparison of numerical and exact solutions is displayed 
in Figure [TTJ We have also calculated the L 1 -norm errors in density for different spatial 
resolutions, where L 1 = J2jAxj\p j — Po(xj)\ with p 0j - and Po(xj) being the numerical and 
exact solutions, respectively. In Table [Q we summarize the results of convergence check for 
this problem. It is thus confirmed that our code has indeed a second order convergence for 
smooth flows (see also Figure [T2|) . 



A. 5. A 2D Isentropic Flow 

We have also performed a 2D computation of the isentropic flow to assess the conver- 
gence rate of our code in a multi-dimensiona l context. The i n itial c ondition for this test 



is the same as in IZhang fc MacFadyenl ( 120061 ); IMorsony et all (120071 ) . The computational 



region is a 2D Cartesian box with 0.0 < x < 3.75 and 0.0 < y < 5.0. The periodic boundary 
condition is adopted for all four sides of the box. The reference state is set to p re f = 1, 
p re f = 100 and v re f = 0. The polytropic EOS with the polytropic constant of K = 100 and 
the adiabatic index of 7 = 5/3 is again employed just as in the ID case. Periodic pulses are 
prepared initially in such a way that they have a spatial period of S = 3.0 in the direction 
given by a unit vector, k = (4/5,3/5), and are uniform in the perpendicular direction. The 
projected spatial periods on the x- and y-directions are 3.75 and 5.0, respectively, which is 
consistent with the size of our computational domain. The initial density profile is given by 
Po(d) (eqs. (lAlj) and (1A2j) ). in which d is a distance from the center of the nearest pulse 
and expressed as d = mod(k • r + 5/2, S) — S/2 with mod(a, b) being a function defined as 
mod(a, b) = a — [a/b\ x b, where [a/b\ denotes the integer part of a/b. The amplitude of the 
pulse is chosen to be a = 1.0 and the width is set to be L = 0.9. The velocity distribution 
in the pulse is determined as in the ID case in the previous subsection so that J_ defined 
for this oblique ID problem should be constant. The simulation is run up to t = 2.4. 

Figure [13] shows the numerical result as a density contour at t — 2.4. In this figure, 
we display the case which the numbers of grid points are set to be 96 and 128 in x— and 
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y— direction, respectively. Just like the one dimensional counter part, the pulse becomes 
asymmetric with the right side of pulse becoming narrower than the left side. We have 
also run some simulations with different numerical resolutions and confirmed again that the 
convergence is approximately of second order (see Tableland Figure IT4"]) . Note that the 
error, Spo, is defined as an L 1 norm in 2D to be 

_ AxAyE jik \p 0jk - p (x j>k )\ 
°Po = a — a v 7 \ \ Ab 

where poj k and po(xj tk ) are the numerical and exact solutions for density at the mesh point 
having an address (j, k). 



A. 6. Dynamical Simulations of Rotational Equilibrium 

In order to check the consistency of the dynamical code with the code of the HSCF 
method employed to construct the rotational equilibrium, we have run long-term simula- 
tions for the stellar envelope in rotational equilibrium, which is obtained by the HSCF 
method (see section I2TT]) . The initial configuration will not change in time if it is indeed in 
dynamical equilibrium. The test will hence validate both the hydrodynamics code with the 
weak gravitational field approximation and the HSCF code simultaneously. The computa- 
tions are essentially the same as those done in Step 3 in the main body except for a different 
treatment of the boundary condition. All the quantities are fixed at the boundaries to the 
values provided by the HSCF calculations unlike in Step 3, where the radial gradients of them 
are artificially reduced to zero. In this subsection, we adopt cgs units. The computational 
domain covers a radial extent of 10 8 cm < r < 2 x 10 10 cm. The simulations are continued 
until t = 100s, which is much longer than the dynamical time scale in the inner region, where 
it is the shortest. Two spatial resolutions have been tried to see the numerical convergence: 
the normal resolution with 230 radial points and higher resolution with 460 points. Just 
as in Step 3, the grid width is determined by the scale height: Ar in = 7.9 x 10 6 cm and 
Ar out = 2.0 x 10 8 cm for the innermost and outermost grids, respectively, for the normal 
resolution, and they are twice finer for the high resolution case. The angular grid is uniform 
and has 60 points in 0° < 9 < 90°. 

Figure [15] shows the density profiles along the rotational axis at the end of the simula- 
tions. The red lines are the profile obtained by the HSCF method and the green ones are 
numerical results. Deviations, which are inevitably induced by the mapping of the initial 
data as well as the difference in the finite-difference methods, are very small and it is indeed 
remarkable that the initial configuration is maintained in such a long time. This is clear 
evidence that both the hydrodynamics code and HSCF code are reliable. It should be also 
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noted that these deviations are even smaller for the high resolution case (see the bottom 
panel of Figure [T5|) . 

A. 7. Sedov- Taylor Problems 

In order to validate our AMR implementation, we have solved the Sedov- Taylor problem. 
Although our code is special relativistic, we use it in the non-relativistic regime here. It is also 
noted that a 2D grid is employed for the computation although the Sedov- Taylor problem is 
one dimensional, since our AMR code used in this paper is specialized for the jet simulation 
and based on the axisymmetric, two-dimensional grid. 

In this simulation, the computational domain, 3 x 10 8 cm < r < 1.8 x 10 10 cm, is covered 
by a uniform mesh with 100 radial grid points (Ar = 1.8 x 10 8 cm) (this section is also 
adopted in cgs units). The internal energy of E = 1.60 x 10 48 erg is deposited initially in the 
central region of r < 6.6 x 10 8 cm. The uniform density is set to be p = lg/cm 3 . We put 
a tiny specific internal energy (e = 10~ 8 x c 2 erg/g) uniformly except for the central region 
mentioned above just for numerical reasons. We adopt the 7-law EOS with 7 = 4/3. We 
impose the free boundary condition both for the inner and outer boundaries. Although all 
9 derivatives vanish initially, they may evolve with time by numerical errors and we set 60 
uniform angular grid points for 0° < 9 < 90°. We impose the axisymmetry and equatorial 
symmetry on the z-axis and equatorial plane, respectively. We vary the resolution of the 
second level mesh to check the numerical convergence. The computations are terminated at 
t = 20s. 

Figure [HI] shows the numerical results (green dots) on the z-axis together with the 
analytical solutions (red lines). From left to right, different resolutions are employed for 
the second level mesh: 3 times, 5 times and 7 times finer in each coordinate direction 
than the first level mesh for the left, middle and right columns, respectively. From top 
to bottom, the profiles of density, pressure and radial velocity on the z-axis are displayed, 
respectively. Note that the horizontal axis is a radius normalized by the radius of the shock 
front, R s h = 1.4 x 10 10 cm. As is evident in this figure, our AMR code successfully reproduces 
the analytical solution with an increasing sharpness of the discontinuity as the second level 
mesh becomes finer. 
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A. 8. A Pulse Propagating through Meshes of Different Refinement Levels 

In the AMR, there is a jump in resolution at the boundary between meshes of different 
refinement levels, which may produce unphysical waves. In order to make sure that the 
effect of the mesh boundary is negligible, we compute an pulse passing through the boundary 
adiabatically. 

The initial pulse profile is the same as that employed in IA.4I Since the present AMR 
code is based on the spherical grid, as mentioned earlier, we need in principle to reformulate 
the problem to accommodate a spherical wave. We can avoid this issue , however, by the 
so-called thin shell approximation, in which the radial range of the computational domain 
is taken to be much smaller than the radius itself and the coordinate curvature can be 
safely ignored. This convenient appr oximatio n has be en w idely used for plane-symmetric 



test problems on the spherical grid lYamadal (see e.g. 119971 ). Here we take r* = 10 and 



Ar* = 1.35 for the representative radius and thickness of the shell, respectively. The initial 
density profile is given by po(r — r*) in Eq. flAip . We assume that all 9 derivatives vanish 
initially. We employ a uniform angular mesh with 10 grid points in 0° < 9 < 0.5°. We 
impose axisymmetry on the z-axis and adopt the free boundary condition at 9 = 0.5°. The 
radial extent of the computational domain is —0.35 + < r < 1 + and is covered by a 
uniform mesh with 210 radial grid points. For the inner region of —0.35 + < r < 0.3 + r*, 
we deploy a second level mesh that is 3 times finer in each coordinate direction than the first 
level mesh. 

Figure [T7I shows the numerical evolution of the pulse. No artificial waves are discernible 
when the pulse passes over the boundary (r = 0.3 + r*) between the meshes of different 
refinement levels (cf. Figure [TT]) . Although the post-pulse state is not identical to the 
reference state, the difference (mainly caused by the grid curvature) is negligible (only ~ 1%). 

In the above run, the pulse is initially located in the region covered by the mesh of higher 
refinement level and moves to the region of lower refinement level. We also run a simulation 
in the opposite case, i.e., the pulse is initially put in the region of lower refinement level 
(or the outer region) and moves to the inner region which is covered by the mesh of higher 
refinement level. This is realized by setting the initial density profile as po(r — r* — 0.7) in 
Eq. ( lAip and determining the velocity distribution so that the Riemann invariant J + = const 
in Eq. ( \A3\i should be constant. The reference state and the amplitude and width of the 
pulse are unchanged. The numerical grids are also the same as above. 

Figure [TBI shows the numerical results for the inward-moving pulse. Just as in the first 
case, the pulse passes through the mesh boundary, producing no discernible artificial wave. 
As expected, the pulse profile at the end of computation is a mirror image of the one for the 
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outward-moving pulse. These test results demonstrate the good behavior of our AMR code 
at the mesh boundary. 

A. 9. Axisymmetric, Relativistic Jet Propagation in a Uniform Matter 

The last test is meant to investigate the effect of AMR resolution on relativistic jet 
propagations and has been frequently used in the literature. The computational domain 
covers the region 0.01 < r < 0.5 and 0° < 9 < 90°. Axisymmetry and equatorial symmetry 
are assumed and a uniform mesh with 150 radial grid points and 60 angular grid points are 
adopted as the first level mesh. We perform two simulations with different AMR resolutions: 
the second level meshes are 3 times and 9 times finer than the first level mesh, respectively. 
The relativistic jet is injected from the inner boundary by hand into the uniform medium. 
The injection parameters are po& = 0.01, pb = 1.70 x 10~ 4 and Vb = 0.99. We employ 7- law 
EOS with 7 = 5/3. The injection parameters give the Mach number of Mb = 6. The half 
opening angle of the jet is chosen as 9h op = 9°. The reference state has a density and pressure 
of Pam = 1-0 and p am = 1.70 x 10~ 4 , respectively. Note that the ambient pressure is the 
same as the jet pressure. The simulation is terminated at t — 2. 

Figure [19] shows the density contours at the end of the computation, t = 2 for the two 
different AMR resolutions. The propagation of the forward shock wave and global struc- 
ture are not very different between the two cases. It is evident that the higher resolution 
captures more complex internal structures. It is well known that the internal jet structure 
never converges as the numerical resolution gets better. This is due to the Kelvin-Helmholtz 
instability, for which the growth rate is greater for shorter wavelengths. Indeed the inter- 
nal structures are very different between the two runs. The obtained numerical results in 
this test are qualitat i vely c onsistent with others in the literature (see e.g. Figure 11 in 



Abdo, A. A., et al. 2009, ApJ, 706, L138 

Abramowicz, M. A., Novikov, I. D., & Paczynski, B. 1991, ApJ, 369, 175 
Aloy, M. A., Miiller, E., Ibanez, J. M., Marti, J. M., 

Blinnikov, S. I., Dunina-Barkovskaya, N. V., & Nadyozhin, D. K. 1996, ApJS, 106, 171 
Campana, S., et al. 2006, Nature, 442, 1008 




REFERENCES 



-23 - 

Del Zanna, L., & Bucciantini, N. 2002, A&A, 390, 1177 
Falk, S. W. 1978, ApJ, 225, L133 
Guiriec, S., et al. 2010, arXiv: 1010.46011 
Hachisu, I. 1986, ApJS, 61, 479 

Ioka, K., Murase, K., Toma, K., Nagataki, S., & Nakamura, T. 2007, ApJ, 670, L77 

Ioka, K. 2010, Progress of Theoretical Physics, 124, 667 

Ito, H., Nagakura, H., Kiuchi, K., Yamada, S., 2011 in preparation 

Kiuchi, K., Nagakura, H., & Yamada, S. 2010, ApJ, 717, 666 

Klein, R. L, & Chevalier, R. A. 1978, ApJ, 223, L109 

Kotake, K., Sato, K., & Takahashi, K. 2006, Reports on Progress in Physics, 69, 971 

Kumar, P., & Narayan, R. 2009, MNRAS, 395, 472 

Kurganov, A., & Tadmor, E. 2000, J. Comput. Phys., 160, 241 

Lazzati, D. 2005, MNRAS, 357, 722 

Lazzati, D., & Begelman, M. C. 2005, ApJ, 629, 903 

Lazzati, D., Morsony, B. J., & Begelman, M. C. 2009, ApJ, 700, L47 

Lazar, A., Nakar, E., & Piran, T. 2009, ApJ, 695, L10 

Lindner, C. C, Milosavljevic, M., Couch, S. M., & Kumar, P. 2010, ApJ, 713, 800 
MacFadyen, A. I., & Woosley, S. E. 1999, ApJ, 524, 262 
MacFadyen, A. L, Woosley, S. E., & Heger, A. 2001, ApJ, 550, 410 
Matzner, C. D., & McKee, C. F. 1999, ApJ, 510, 379 

Milosavljevic, M., Lindner, C. C, Shen, R., & Kumar, P. 2010. larXiv: 1007.07631 
Mizuta, A., Yamasaki, T., Nagataki, S., & Mineshige, S. 2006, ApJ, 651, 960 
Mizuta, A., & Aloy, M. A. 2009, ApJ, 699, 1261 
Mizuta, A., Nagataki, S., & Aoi, J. 2010. larXiv:1006.2440l 



-24- 



Mizuta, A., Kino, M., & Nagakura, H. 2010, ApJ, 709, L83 

Morsony, B. J., Lazzati, D., & Begelman, M. C. 2007, ApJ, 665, 569 

Nagakura, H., & Yamada, S. 2008, ApJ, 689, 391 

Nagakura, H., Itoh, H., Kiuchi, K., Yamada, S., 2011 in preparation 

Paczynski, B. 1998, ApJ, 494, L45 

Pe'er, A. 2008, ApJ, 682, 463 

Piran, T. 2004, Reviews of Modern Physics, 76, 1143 

Pons, J. A., Ma Marti, J., M&uumlller, E. 2000, Journal of Fluid Mechanics, 422, 125 

Rezzolla, L., & Zanotti, O. 2002, Physical Review Letters, 89, 114501 

Ryde, F., et al. 2010, ApJ, 709, L172 

Soderberg, A. M., et al. 2008, Nature, 454, 246 

Toma, K., Wu, X.-F., & Meszaros, P. 2010. larXiv:1002.2634l 

Tominaga, N., Maeda, K., Umeda, H., Nomoto, K., Tanaka, M., Iwamoto, N., Suzuki, T., 

Waxman, E., Meszaros, P., & Campana, S. 2007, ApJ, 667, 351 

Woosley, S. E. 1993, ApJ, 405, 273 

Woosley, S. E., & Bloom, J. S. 2006, ARA&A, 44, 507 

Woosley, S. E., & Heger, A. 2006, ApJ, 637, 914 

Yamada, S. 1997, ApJ, 475, 720 

Yonetoku, D., Murakami, T., Nakamura, T., Yamazaki, R., Inoue, A. K., & Ioka, K. 2004, 
ApJ, 609, 935 

Zhang, W., Woosley, S. E., & MacFadyen, A. I. 2003, ApJ, 586, 356 
Zhang, W., Woosley, S. E., & Heger, A. 2004, ApJ, 608, 365 
Zhang, W., & MacFadyen, A. I. 2006, ApJS, 164, 255 



This preprint was prepared with the A AS IATgX macros v5.2. 



Table 1. Numerical Errors for Different Resolutions in the ID Isentropic-Flow Problem 



Number of Grid Points 


L 1 Error (%) 


Convergence Rate (a) 


100 


0.58 




200 


0.17 


1.83 


400 


3.67E-2 


2.17 


800 


7.96E-3 


2.20 


1600 


1.80E-3 


2.14 


3200 


4.11E-4 


2.13 


6400 


9.60E-5 


2.10 



Note. — The errors of density are evaluated at t = 0.8. In the 
right-most column, the powers in the expression L 1 oc N~ a , where 
N denotes the number of grids, are given. 
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Table 2. Numerical Errors for Different Resolutions in the 2D Isentropic-Flow Problem 



Number of Grid Points 


Spo (%) 


Convergence Rate (a) 


48 x 64 


0.90 




96 x 128 


0.24 


1.93 


192 x 256 


6.17E-2 


1.93 


384 x 512 


1.24E-2 


2.32 


768 x 1024 


2.70E-3 


2.20 



Note. — The errors of density are evaluated at t = 2.4. See 
the text for the definition of Spg. In the right-most column, 
the powers in the expression Spo oc N~ a , where N denotes the 
number of grids, are given. 
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Fig. 1. — The density profiles on the rotational axis (z-axis) and equato r (x-axis) for the 
envelope model in this paper and model 16TI by IWoosley fc Hegerl (120061 ) (left panel) and 
the density contour (log scale) in the meridian section of the same envelope model (right 
panel) . 
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Fig. 2. — The profile of specific angular momentum on the equatorial plane for the envelope 
model in this paper as a function of included mass (green line). The red line shows the 
specific angular momenta at the innermost stable circular orbit (ISCO) for Schwartzschild 
Black Holes function of their masses. 
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Fig. 3. — The profiles of density (upper), pressure (middle), and radial velocity as well 
as sound velocity (bottom) for the wind model employed in this paper. The red lines in 
the upper and middle panels show power-laws for comparison. As shown in bottom panel, 
although the radial velocity is initially subsonic, it passes a sonic surface, which is located 
at r ~ 7.5 x 10 10 cm. 
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Fig. 4. — The time evolution of the density profile on the equator during the envelope 
collapse (top panel) and the density contours in the meridian section at the time of the 
breakout of the shock wave produced by centrifugal bounce (bottom left panel) and ~ 10s 
later (bottom right panel). 
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Fig. 5. — The density contours (upper four panels) and the time evolutions of the Lorentz 
factors on the rotational axis (lower four panels) for model M20s (left column) and model 
M50s (right column). The time, t, is measured from the instant, at which the jet is injected. 
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Fig. 6. — The light curves of photospheric emissions (top panels), the evolutions of the 
observed temperature (middle panels) and of the photospheric radius (bottom panels) as a 
function of the observed time for model M20s (left column) and model M50s (right column). 
See the body for details. 
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Fig. 7. — Numerical results (dots) for the ID relativistic shock tube problems without 
tangential velocities. The rest mass density (p), pressure p and velocity v x are shown. 
The exact solutions (solid lines) are also displayed for comparison. The left (right) panel 
corresponds to case 1 (case 2) at t = 0.4 (t = 0.35). 
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Fig. 8. — Numerical results (dots) for the ID relativistic shock tube problems with tangential 
velocities. A uniform mesh with 400 grid points are employed. The exact solutions (solid 
lines) are also displayed for comparison. We change Vy from left to right as Vy = 0, 0.9, 0.99 
and from top to bottom as Vy = 0, 0.9, 0.99. The den sity, pressure and x-component of 
velocity are shown. See Figure 15 in iMizuta et al.l ( 120061 ) for comparison. 
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Fig. 9. — Numerical results (dots) for the same problem as in the previous figure (Vy, Vy) = 
(0.9, 0.9) with different resolutions. These are meant to check the numerical convergence. 
The exact solutions (solid lines) are displayed for comparison. The left panels show the 
rest mass density and pressure, whereas the middle (right) panels display the x-component 
(y-component )of velocity. From top to bottom, the number of grid points are 800, 1600, 
3200 and 6400, respectively. 
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Fig. 10. — The contour plot of rest mass density in the logarithmic scale at t — 0.4 for the 
2D Riemann problem. 
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Fig. 11. — The initial (t = 0) and simulated (t = 0.8) density (top), pressure (middle), 
velocity (bottom) profiles for the ID isentropic flow together with the exact solution (solid 
lines). A uniform mesh with 400 grid points is employed. 
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Fig. 12. — The L 1 density errors for the ID isentropic flow as a function of the number of 
grid points. 
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Fig. 13. — The numerical result of the density structure for the 2D isentropic-flow problem. 
The spatial resolution is effectively equivalent to the numbers of grid points of (96,128) in 
the x— and y— directions. 
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Fig. 14. — The density errors Spo for the 2D isentropic flow as a function of the number of 
grid points. 
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Fig. 15. — The density profiles along the rotational axis at t — 100s after the long-term 
dynamical simulations of the envelope in rotational equilibrium. Two spatial resolutions are 
employed with the upper panel showing the result for 230 radial grid points (green line) 
together with the initial profile (red line), whereas the middle panel corresponding to the 
result for 460 points. 




Fig. 16. — The computed profiles of density (top), pressure (middle) and radial velocity 
(bottom) together with the exact solution (red lines) for the Sedov- Taylor problem. The 
AMR resolution becomes higher from left to right. 
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Fig. 17. — The computed profiles of density (top), pressure (middle) and radial velocity 
(bottom) at t = 0, 0.4, 0.8 for the right-moving pulse. The mesh boundary is located at 
r - 10 4 = 0.3. 
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Fig. 18- 



The same as Figure [T7] but for the left-moving pulse. 
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Fig. 19. — Axisymmetric, relativistic jet propagations in a uniform medium. The density 
contours at t = 2 are displayed. The left panel shows the result for the case, in which the 
second level mesh is 3 times finer than the first level mesh and the right panel gives the 
result when a 9 times finer second level mesh is employed. 



